      SUBROUTINE HALLEY(TIME1,POS)
C
C SUBROUTINE HALLEY TAKES A TIME AND RETURNS THE POSITION OF HALLEY'S 
C COMET WITH RESPECT TO EARTH CENTERED COORDINATES IN MEAN EQUATOR AND 
C EQUINOX OF 1950.0
C
C*****************************************************************
C
C THE METHOD FOR FINDING THE POSITION OF HALLEY'S COMET IN THIS
C SUBROUTINE IS AS FOLLOWS: THE SIX KEPLERIAN ORBITAL ELEMENTS
C OF HALLEY'S COMET ARE DEFINED IN THE DATA STATEMENT.
C
C*****************************************************************
C
C THIS SUBROUTINE WAS CODED BY BARBARA H. ROBERTS, 6/81.
C UPDATED............
C            1/83. C PETRUZZO. COMMENT CHANGES, REPLACE ELREC WITH 
C                  TOCART, AND ASSOCIATED CHANGES.
C            5/85. S LAMBROS. UPDATE ORBITAL ELEMENTS AND TIMEP, AS
C                  PER D.K. YEOMANS OF JPL. ALSO CHANGED ASSOCIATED
C                  COMMENTS.
C******************************************************************
C
C VARIABLE DESCRIPTION  DIMENSION     I/O    TYPE     DESCRIPTION
C
C        TIME1              1          I      R*8     SEC SINCE 1/1/50,
C                                                     0.0 HOURS
C
C        POS                3          O      R*8     POSITION OF
C                                                     COMET IN X,Y,Z
C                                                     COORDINATES WRT 
C                                                     EARTH IN KM
C 
C INTERNAL VARIABLES:
C
C        ELEM               6         N/A     R*8     KEPLERIAN ELEMENTS
C                                                     (AU,DEG) OF HALLEY
C
C        U                  1         N/A     R*8     GRAVITATIONAL
C                                                     CONSTANT
C
C        ELEM(6)            1         N/A     R*8     MEAN ANOMALY
C        
C**********************************************************************
C
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 AU/149.6D6/
      DIMENSION ELEM(6),POS(3),COEF(4),RHALEY(3),VHALEY(3),SUNVEC(3)
C
C   HALLEY ORBITAL ELEMENTS. THESE ARE FROM D. K. YEOMANS OF JPL -
C   HIS LATEST CALCULATIONS AS OF 5/17/85. UNITS ARE A.U. AND DEG.,
C   AND ARE REFERRED TO THE ECLIPTIC AND EQUINOX OF 1950.0.
C   PERIHELION TIME(=TIMEP) CORRESPONDS TO 1986, FEB. 9.45069 E.T.
      DATA ELEM/17.9390877D0, 0.9672726D0,
     1  162.23931D0, 58.14398D0, 111.84677D0, 0.D0/,   ! DEGREES
     2  TIMEP/2446470.95069D0/
C
C   COEFFICIENTS FOR OBLIQUITY COMPUTATION.
      DATA COEF/23.4522944D0, -.0130125D0,-.1638889D-5,-.5027778D-6/
C
      REAL*8 JD50/2433282.5D0/
      REAL*8 DEGRAD/ 57.29577951308232D0 /
      REAL*8 TWOPI/ 6.283185307179586D0 /
      INTEGER INIT/0/
C 
C
      IF(INIT.NE.0) GO TO 55
      INIT=1
C    OBLIQUITY AT 1/1/50,  0.0 HRS.
      TT=0.5D0
      TILT=(COEF(1)+COEF(2)*TT+COEF(3)*TT**2+COEF(4)*TT**3)/DEGRAD
      SINI=DSIN(TILT)
      COSI=DCOS(TILT)
C    SUN MU IN KM**3/SEC**2
      U=1.32715445D11
C    SUN MU IN AU**3/SEC**2
      U=U/AU**3
C    PERIOD IN DAYS.
      PERIOD=TWOPI*DSQRT(ELEM(1)**3/(U*86400.D0**2))
C    TPERIH=PERIHELION TIME IN DAYS SINCE 1/1/50, 0.0 HRS E.T.
      TPERIH= TIMEP-JD50
C    HALLEY ANGLE ELEMENTS CONVERTED TO RADIANS.
      DO 54 I=3,6
   54 ELEM(I)=ELEM(I)/DEGRAD
   55 CONTINUE
C
C    FIND RHALEY, THE POSITION VECTOR OF HALLEY WRT THE SUN AND IN 
C    ECLIPTIC COORDINATES, MEAN OF 1950.0.
      DT=TIME1/86400.0D0-TPERIH
      ELEM(6)=DT/PERIOD * TWOPI
      CALL TOCART(U,ELEM,0,RHALEY,VHALEY,6,IERR)
C
C    NOW ROTATE RHALEY TO EARTH'S EQUATOR
      POS2=RHALEY(2)*COSI-RHALEY(3)*SINI
      POS3=RHALEY(3)*COSI+RHALEY(2)*SINI
      RHALEY(2)=POS2
      RHALEY(3)=POS3
C
C    FIND THE POSITION OF SUN IN MEAN OF 1950.0.
      CALL SOLM50(TIME1,SUNVEC,2)
C
C    NOW GET AND RETURN HALLEY POSITION IN GEOCENTRIC, MEAN EQUATOR
C    AND EQUINOX OF 1950.0.  UNITS ARE KM.
      DO 300 I=1,3
  300 POS(I)=RHALEY(I)*AU+SUNVEC(I)
C
      RETURN
      END
